#include <cmath>
#include <omp.h>

#include "config.h"
#include "network.h"
#include "math_op.h"

#include "nef.h"

size_t axis_length, ne_nf_length;
size_t single_axis_length, single_ne_nf_length;

double *axis_c, *single_axis_c; // C stands for concentration.
size_t *axis_base;
double *NE_NF_results, *single_NE_NF_results;
std::deque<int> *ec_index_lists;
double **single_ec_lists;

#define ID_MAX   (0)
#define ID_MIN   (1)
#define ID_DELTA (1)
static double NE_NF_Min_Max_results[2][2];
static int **single_ec_index_lists;

void ne_nf_prepare_space()
{
    //Calculate the axis length.
    axis_length = static_cast<size_t>((c_end - c_start) / c_step + 1.0);
    //Based on axis length, calculate the nef result space.
    ne_nf_length = static_cast<size_t>(pow(axis_length, num_of_compounds));
    //Calculate the concentration axis mapping index.
    axis_c = new double[axis_length];
    #pragma omp parallel for
    for(int i=0; i<axis_length; ++i)
    {
        // Calculate the axis concentration.
        axis_c[i] = c_start + c_step * i;
    }
    //Allocate the NE/NF result array.
    NE_NF_results = static_cast<double *>(
                malloc(sizeof(double) * ne_nf_length * 2));
    //Calculate the division for each axis.
    axis_base = new size_t[num_of_compounds];
    for(size_t i=0; i<num_of_compounds; ++i)
    {
        //Calculate the axis base.
        axis_base[i] = static_cast<size_t>(pow(axis_length,
                                               num_of_compounds - i - 1));
    }
    //Create deque for all the EC list.
    ec_index_lists = new std::deque<int>[ec_list.size()];

}

void ne_nf_prepare_axis_space()
{
    //Calculate the single compound length.
    single_axis_length = static_cast<size_t>((s_end - s_start) / s_step + 1.0);
    //Based on the single axis_length, calculate the space it needs.
    single_ne_nf_length = single_axis_length * num_of_compounds;
    //Prepare the axis.
    single_axis_c = new double[single_axis_length];
    //Allocate the NE/NF result array.
    single_NE_NF_results =
            static_cast<double *>(malloc(sizeof(double) * single_ne_nf_length * 2));
    //Calculate the axis concentrations.
    #pragma omp parallel for
    for(int i=0; i<single_axis_length; ++i)
    {
        //Calculate the position concentration.
        single_axis_c[i] = s_start + s_step * i;
    }
    //Prepare the EC index list.
    single_ec_index_lists = new int *[num_of_compounds];
}

void global_index_to_concentration(size_t residual, double *cc_coordinate)
{
    //According to the index i, calculate the space coordinate.
    size_t axis_val;
    size_t *coordinate = new size_t[num_of_compounds];
    for(size_t j=0; j<num_of_compounds; ++j)
    {
        //Calculate each axis.
        axis_val = residual / axis_base[j];
        //Save the axis val.
        coordinate[j] = axis_val;
        //Remove the axis from residual.
        residual -= axis_val * axis_base[j];
    }
    // Based on the axis, read the value.
    for(size_t j=0; j<num_of_compounds; ++j)
    {
        //Based on the axis concentration, read the data.
        cc_coordinate[j] = axis_c[coordinate[j]];
    }
    //Clear the coordinate memory.
    delete[] coordinate;
    // Check the global log switch.
    if(c_log)
    {
        // Changing all the coordinator into power value (reverse of log).
        for(size_t j=0; j<num_of_compounds; ++j)
        {
            //Based on the axis log concentration,
            //calculate the power result.
            cc_coordinate[j] = pow(10, cc_coordinate[j]);
        }
    }
}

inline void calculate_ev_by_concentration(double *ev, double *cc_coordinate)
{
    //Convert cc coordinate unit (from Mol to uMol).
    for(size_t i=0; i<num_of_compounds; ++i)
    {
        cc_coordinate[i] /= 1000000.0;
    }
    //Loop for all the nodes.
    for(size_t i=0; i<num_of_nodes; ++i)
    {
        //Based on the value of Ki, calculate the ev.
        //Check the first value of the node.
        if(compound_Ki_list[0][i] < 0.0)
        {
            //This node is ignored, set the ev as CONNECTED constant.
            ev[i] = EV_CONNECT;
            continue;
        }
        //Sum up all compounds Ki.
        double base = 1.0;
        for(size_t j=0; j<num_of_compounds; ++j)
        {
            base += cc_coordinate[j] / compound_Ki_list[j][i];
        }
        //Set the EV value.
        ev[i] = pow(10, 2.303 * (1-1/base));
    }
}

inline void calculate_single_ev_by_concentration(double *ev, double cc, int compound_index)
{
    //Convert cc into correct unit (from Mol to uMol).
    cc /= 1000000.0;
    //Get the index.
    const std::vector<double> &Ki_list = compound_Ki_list[compound_index];
    //Loop for all the nodes.
    for(size_t i=0; i<num_of_nodes; ++i)
    {
        //Check the target compound index value.
        if(Ki_list.at(i) < 0.0)
        {
            //This node is ignored, set the ev as CONNECTED constant.
            ev[i] = EV_CONNECT;
            continue;
        }
        //Calculate the target ev.
        ev[i] = pow(10, 2.303 * (1-1/(1.0+cc/Ki_list.at(i))));
    }
}

//Priority queue structure.
typedef struct QueueNode
{
    int vid;
    double value;
    QueueNode(int mvid=0, double mvalue=0.0) :
        vid(mvid),
        value(mvalue)
    {}

    bool operator < (const QueueNode& rhs) const {
      return value > rhs.value;
    }
} QueueNode;

void search_ne_nf(double *ev, double *ne_nf)
{
    //Now we have the EV value, based on the network value, find NE and NF.
    //Opt Dijkstra.
    double NE = 0.0, NF = 0.0;
    bool *final = new bool[num_of_nodes];
    double *dist = new double[num_of_nodes];
    QueueNode current;
    for(int m=0; m<num_of_nodes; ++m)
    {
        // Execute Dijkstra for all vertex.
        for(int i=0; i<num_of_nodes; ++i)
        {
            final[i] = false;
            dist[i] = (i == m) ? 0.0 : EV_INFINITY;
        }
        // Prepare the priority queue.
        std::priority_queue<QueueNode> queue;
        // Put the current node to queue.
        current.vid = m;
        current.value = 0.0f;
        queue.push(current);
        // Loop until all the queue is executed.
        while(!queue.empty())
        {
            // Pop the top from the queue.
            current = queue.top();
            queue.pop();
            // Check final, although I think it is unnecessary.
            if(final[current.vid])
            {
                // Do not need to search again.
                continue;
            }
            // Mark as final.
            final[current.vid] = true;
            // Loop for checking all the edges.
            for(int i=0; i<edges[current.vid].size(); ++i)
            {
                //Get the target id.
                int to_id = edges[current.vid][i];
                // Check final of the vid point.
                if(!final[to_id])
                {
                    //Try the value.
                    double trial = ev[current.vid] + dist[current.vid];
                    if(dist[to_id] > trial)
                    {
                        // Update the distance.
                        dist[to_id] = trial;
                        // Add the node to queue.
                        queue.push(QueueNode(to_id, trial));
                    }
                }
            }
        }
        // Calculate NE and NF.
        // NE
        for(int i=0; i<num_of_nodes; ++i)
        {
            if(dist[i]!=0 && !is_close(dist[i], EV_INFINITY))
            {
                NE += 1.0 / dist[i];
            }
        }
        // NF.
        for(auto i : vertex_exit)
        {
            //Add the exit node to NF, only when m is not i.
            if(m != i)
            {
                if(dist[i]!=0 && !is_close(dist[i], EV_INFINITY))
                {
                    NF += 1.0 / dist[i];
                }
            }
        }
    }
    delete[] final;
    delete[] dist;
    //Save NE and NF to result.
    ne_nf[0] = NE;
    ne_nf[1] = NF;
}

void ne_nf_calculate()
{
    //For all the concentration, calculate its NE and NF.
    // Disable the macro to enable serial processing.
    #pragma omp parallel for
    for(int i=0; i<ne_nf_length; ++i)
    {
        //Prepare the compound concentration coordinate.
        double *cc_coordinate = new double[num_of_compounds];
        //Calculate the concentration coordinate.
        global_index_to_concentration(i, cc_coordinate);
        //Based on the coordination, calculate the network edge weight.
        double *ev = new double[num_of_nodes];
        calculate_ev_by_concentration(ev, cc_coordinate);
        delete[] cc_coordinate;
        //Based on EV, calculate the NE and NF.
        search_ne_nf(ev, NE_NF_results + i * 2);
        delete[] ev;
    }
}

void ne_nf_single_calculate()
{
    //Loop for all the axis value.
    for(int m=0; m<num_of_compounds; ++m)
    {
        //Get the destination address.
        double *single_result = single_NE_NF_results + m * single_axis_length * 2;
        //For each compounds, we are going to calculate its NE/NF.
        //For all the concentration, calculate its NE and NF.
        // Disable the macro to enable serial processing.
        #pragma omp parallel for
        for(int i=0; i<single_axis_length; ++i)
        {
            //Based on the concentration on axis, calculate EV.
            double *ev = new double[num_of_nodes];
            calculate_single_ev_by_concentration(ev, single_axis_c[i], m);
            //Based on EV, calculate the NE and NF.
            search_ne_nf(ev, single_result + i * 2);
            delete[] ev;
        }
    }
}

void ne_nf_min_max_calculate()
{
    //Prepare the min and max edge value.
    double *buff_ev_min_max = new double[num_of_nodes * 2];
    //Prepare the pointer value.
    double *ev_min_max[2];
    ev_min_max[ID_MAX] = buff_ev_min_max + num_of_nodes * ID_MAX;
    ev_min_max[ID_MIN] = buff_ev_min_max + num_of_nodes * ID_MIN;
    //Assign edge value to buff min and max.
    for(size_t i=0; i<num_of_nodes; ++i)
    {
        //NEmax and NEmin are the maximal and minimal NE when all EVs are set as
        //1 and 200 -- p.8
        ev_min_max[ID_MAX][i] = EV_CONNECT;
        ev_min_max[ID_MIN][i] = EV_MAX_RESIST;
    }
    //Calculate the result of min and max.
    // Disable the macro to enable serial processing.
    #pragma omp parallel for
    for(int i=0; i<2; ++i)
    {
        //Execute the search function.
        search_ne_nf(ev_min_max[i], NE_NF_Min_Max_results[i]);
    }
    //Clear the buffer.
    delete[] buff_ev_min_max;
}

void ne_nf_to_nef()
{
    //Calculate the NEmax-NEmin, NFmax-NFmin
    const double NEdelta = NE_NF_Min_Max_results[ID_MAX][0] - NE_NF_Min_Max_results[ID_MIN][0],
                 NFdelta = NE_NF_Min_Max_results[ID_MAX][1] - NE_NF_Min_Max_results[ID_MIN][1],
                 NEmax = NE_NF_Min_Max_results[ID_MAX][0],
                 NFmax = NE_NF_Min_Max_results[ID_MAX][1];
    //Loop and convert all the NE and NF to NEF.
    #pragma omp parallel for
    for(int i=0; i<ne_nf_length; ++i)
    {
        //Get the NE/NF pointer.
        double *ne_nf_block = NE_NF_results + i * 2;
        //Calculate the NEF value.
        double NEF = sqrt(((NEmax - ne_nf_block[0]) / NEdelta) * // NEd
                          ((NFmax - ne_nf_block[1]) / NFdelta)   // NFd
                );
        //Save the NEF value to NE position.
        ne_nf_block[0] = NEF;
    }
}

void ne_nf_single_to_nef()
{
    //Calculate the NEmax-NEmin, NFmax-NFmin
    const double NEdelta = NE_NF_Min_Max_results[ID_MAX][0] - NE_NF_Min_Max_results[ID_MIN][0],
                 NFdelta = NE_NF_Min_Max_results[ID_MAX][1] - NE_NF_Min_Max_results[ID_MIN][1],
                 NEmax = NE_NF_Min_Max_results[ID_MAX][0],
                 NFmax = NE_NF_Min_Max_results[ID_MAX][1];
    //Loop for all the compounds, convert their NE/NF value.
    for(int m=0; m<num_of_compounds; ++m)
    {
        //Get the destination address.
        double *single_result = single_NE_NF_results + m * single_axis_length * 2;
        //Loop for all the result.
        #pragma omp parallel for
        for(int i=0; i<single_axis_length; ++i)
        {
            //Get the NE/NF pointer.
            double *ne_nf_block = single_result + i * 2;
            //Calculate the NEF value.
            double NEF = sqrt(((NEmax - ne_nf_block[0]) / NEdelta) * // NEd
                              ((NFmax - ne_nf_block[1]) / NFdelta)   // NFd
                    );
            //Save the NEF value to NE position.
            ne_nf_block[0] = NEF;
        }
    }
}

void filter_ec_from_nef()
{
    //Loop and search for all the NEF data.
    #pragma omp parallel for
    for(int m=0; m<ec_list.size(); ++m)
    {
        //Buffer the value to compare.
        double target_ec = ec_list.at(m);
        //Based on the ec value, check search in the NEF result.
        for(int i=0; i<ne_nf_length; ++i)
        {
            //Get the NEF pointer.
            if(is_close(target_ec, (NE_NF_results + i * 2)[0], ec_epsinon))
            {
                //Append the NEF index result to queue.
                ec_index_lists[m].push_back(i);
            }
        }
    }
}

void single_ec_index_to_concentration()
{
    //Allocate all the concentration array.
    single_ec_lists = new double *[num_of_compounds];
    //Convert all the index of the array to concentration.
    for(int ci=0; ci<num_of_compounds; ++ci)
    {
        //Prepare the memory.
        single_ec_lists[ci] = new double[ec_list.size()];
        //For all the concentration in the list, loop to convert all.
        #pragma omp parallel for
        for(int m=0; m<ec_list.size(); ++m)
        {
            if(single_ec_index_lists[ci][m] < 0)
            {
                //Set a constant minus number to that position.
                single_ec_lists[ci][m] = -1.0;
            }
            else
            {
                //Look up the table.
                single_ec_lists[ci][m] = single_axis_c[single_ec_index_lists[ci][m]];
            }
        }
        //Recovert the ec index array.
        delete[] single_ec_index_lists[ci];
    }
    //Recover the ec index array.
    delete[] single_ec_index_lists;
}

void filter_ec_single_from_nef()
{
    bool error_message = false;
    //Loop and search for all the NEF data for each compounds.
    for(int ci=0; ci<num_of_compounds; ++ci)
    {
        //Calculate the compound NEF destination address.
        double *single_result = single_NE_NF_results + ci * single_axis_length * 2;
        //Prepare the vector.
        single_ec_index_lists[ci] = new int[ec_list.size()];
        //Loop and search for all the NEF data.
        #pragma omp parallel for
        for(int m=0; m<ec_list.size(); ++m)
        {
            //Buffer the value to compare.
            double target_ec = ec_list.at(m);
            //Prepare the EC data list.
            std::deque<std::pair<int, double>> compound_ec;
            //Based on the EC value, get all the possible NEF data here.
            for(int i=0; i<single_axis_length; ++i)
            {
                //Get the NEF pointer.
                double nef_value = (single_result + i * 2)[0];
                if(is_close(target_ec, nef_value, ec_epsinon))
                {
                    //Append the NEF index result to queue.
                    compound_ec.push_back(std::pair<int, double>(i, nef_value));
                }
            }
            //Check the size of compound EC.
            if(compound_ec.size()==0)
            {
                error_message=true;
                printf("\nEC%d is failed to be found in compound No.%d",
                       static_cast<int>(target_ec*100),
                       ci);
                //Set the position with a -1.
                single_ec_index_lists[ci][m] = -1;
                //Skip.
                continue;
            }
            //Calculate the middle value of the deque.
            int middle_pos = compound_ec.size() / 2;
            std::nth_element(
                        compound_ec.begin(), compound_ec.begin()+middle_pos,
                        compound_ec.end(),
                        [](std::pair<int, double> left, std::pair<int, double> right)
            {
                return left.second < right.second;
            });
            //Now, the target position is the correct value, get it!
            single_ec_index_lists[ci][m] = compound_ec.at(middle_pos).first;
        }
    }
    if(error_message)
    {
        printf("\n");
    }
    //Convert all the ec index to concentration.
    single_ec_index_to_concentration();
}
